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ABSTRACT 

Aims. An alternative to the traditional method for modeling kinematics of the Earth's rotation is proposed. The purpose of developing 
the new approach is to provide a self-consistent and simple description of the Earth's rotation in a way that can be estimated directly 
from observations without using intermediate quantities. 

Methods. Instead of estimating the time series of pole coordinates, the UTl-TAI angles, their rates, and the daily offsets of nutation, 
it is proposed to estimate coelficients of the expansion of a small perturbational rotation vector into basis functions. The resulting 
transformation from the terrestrial coordinate system to the celestial coordinate system is formulated as a product of an a priori 
matrix of a finite rotation and an empirical vector of a residual perturbational rotation. In the framework of this approach, the specific 
choice of the a priori matrix is irrelevant, provided the angles of the residual rotation are small enough to neglect their squares. The 
coefficients of the expansion into the B-spline and Fourier bases, together with estimates of other nuisance parameters, are evaluated 
directly from observations of time delay or time range in a single least square solution. 

Results. This approach was successfully implemented in a computer program for processing VLBI observations. The dataset from 
1984 through 2006 was analyzed. The new procedure adequately represents the Earth's rotation, including slowly varying changes in 
UTl-TAI and polar motion, the forced nutations, the free core nutation, and the high frequency variations of polar motion and UTl. 
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1. Introduction 

, As we learn from an elementary physics course in high school, 

■ the rotation of a rigid body can be represented by three angles. 
' Euler angles are usually selected as parameters, although it is not 

the only choice. We can say that rotation of an arbitrary body is 
defined if a functional dependence of Euler angles on time is 
known. 

However, when we are dealing with the Earth, the tradi- 
tional way of representing the kinematics of the Earth's rota- 
, tion is not as simple. In modern textbooks that follow t he for - 

■ malism of Newcomb-Andoyer, for instance. ISeidelmannI (1 1 992h . 

■ the Earth's rotation is represented as a product of 9 matrices 
that depend on many parameters. Analytical expressions are 
provided for some of them, but other parameters are supposed 
to be determined from observations. These parameters are de- 
fined though intermediate quantities. These intermediate quan- 
tities like the Earth's rotation axis, the celestial ephemeris pole, 
the true equinox of date, the non-rotating origin and others, are 
not objects that can be observed, but idealistic concepts that do 
not have a clear, intuitive interpretation. 

According to the traditional approach, (e.g. 
[McCarthy & Petit, eds] 12004*). in order to get the Earth's 
orientation at any given moment, one should first compute 
the values of intermediate parameters that have more than a 
thousand terms, then interpolate tables of corrections produced 
by smoothing results of analysis of observations, and finally 
compute the product of rotational matrices that depend on these 
intermediate quantities. 

At first glance it seems that the complexity of represent- 
ing the Earth's rotation is unavoidable, since it should reflect 
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the complexity of a theory for modeling the Earth's rotation 
at a level of the accuracy of observations, i.e. on the order of 
10""* racQ. This would be true if theoretical models were pre- 
cise to that level. Roughly speaking, the Earth's rotation can be 
considered as consisting of two components, the tidally driven 
component with precisely known frequencies and the compo- 
nent driven by an exchange of the angular momentum between 
the solid Earth and geophysical fluids. The latter component is 
not predictable in principle. The atmosphere contributes to the 
UTl at a level of 10"^ rad, more than three orders of magnitude 
higher than the accuracy of observations. For a long period of 
time it was considered that the quasi-diurnal motion of the pole, 
namely the precession and nutation, can be modeled with a pre- 
cision comparable to the accura cy of observat i ons. F irst results 
of VLBI analysis presented by iHerring et alj ( Il986l) shattered 
this belief. It was discovered that even the most advanced nuta- 
tion theory of Wahr ( 1980) was not accurate enough. Numerous 
attempts to build a theory of tidally driven quasi-diurnal mo- 
tion were made, but they were not completely successful. In 
order to reduce the disagreement between theories and obser- 
vations, the authors had to resort to adjusting some parameters 
of their theories to quantities deri ved from the same observa- 
tions of the quasi-diurnal motion jMathews et al.l ll99lL 120021 : 
iGetino & Ferrandizll200l l). It was also soon realized that there 
are two other constituents of the quasi-diurna l polar m otion at a 
level of 1 nanoradian: the free core nutation (Moritz 1987) and 
the nutation excited by the atmosphere (Bizouard et al. 199^. 
These constituents currently cannot be predicted and, presum- 
ably, they cannot be predicted in principle. Therefore, even if a 
precise theory of forced nutation is built in the future, one should 

' SI units are used throughout the paper. Conversion factors to non- 
standard units: I • 10^'rad » 0.2Imas a; 14//sec 
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apply parameters determined from observations in order to rep- 
resent the quasi-diurnal motion. 

Recognizing that both components in the Earth rotation can- 
not be predicted with accuracy comparable to the precision of 
observations, prompts us to reconsider approaches to represent- 
ing the Earth's rotation and the role of the theory. The quasi- 
diurnal motion should be described with the use of parameters 
determined by continuous observations in a similar way as the 
UTl and Chandler polar motion. The theory of nutation should 
be considered not as a tool for data reduction or for predicting 
the Earth's orientation, but as means for validating geophysical 
models. At the same time a theory of the Earth's rotation can 
provide valuable guidance for building empirical mathematical 
models. 

The goals of this paper is to build such an empirical mathe- 
matical model, to demonstrate using a long dataset of observa- 
tions that it is feasible and to show that this approach describes 
the Earth rotation at least as well as the traditional way. It should 
be noted that an empirical mathematical model has a different 
meaning than a theoretical model. The theoretical model relates 
a function of time that describes the Earth rotation with specific 
properties of the Earth's body in the form of a solution for the 
equation of dynamics. The empirical mathematical model relates 
this function of time to observations using a parameter estima- 
tion technique. A minimal requirement for an empirical model 
is to represent the phenomena with the least possible errors for 
the entire interval of observations and to provide estimates of 
uncertainties. We will also try to satisfy two additional require- 
ments: the model should be simple and the parameters of the 
model should not be strongly correlated. If parameter estimates 
are strongly correlated, their comparison with theoretical predic- 
tions becomes problematic. We will represent the Earth orienta- 
tion parameters (EOP) in the form of an expansion over a family 
of basis functions. 

The procedure for developing an empirical model of the 
Earth's rotation is presented in the rest of the paper The choice 
of the a priori model and basis functions is described in section|2] 
the proposed mathematical models is described in Sect. [3] the 
strategy of analysis of the 22 year long dataset of VLBI obser- 
vations is presented in Sect. |4l and the results of solutions are 
discussed in Sect.|5] Concluding remarks are given in Sect.|6] 

2. Choice of the a priori model and the basis 
functions 

2.1. Model of observations 

We consider here that stations observe K celestial physical 
bodies. It is assumed that each station is associated with a ref- 
erence point. In the case of VLBI antennas with intersecting 
axes, this is the intersection point of the axes. Observing sta- 
tions receive electromagnetic radiation emitted by celestial bod- 
ies, and each sample of the received signal is associated with a 
time stamp from a local frequency standard synchronized with 
the GPS time. Analysis of voltage and time stamps of received 
radiation eventually allows us to derive photon propagation time 
from reference points of observed bodies to reference points of 
observing stations. These distances depend on the relative posi- 
tions of stations with respect to the observed bodies. The instan- 
taneous coordinate vector of station ;, r,(f), at a given moment 
of time is represented as the sum of a rotation and translation 
applied to a vector r,(r()) at initial epoch to as 

r,(f) = VH(f)r,(f(,) + r(f) + rf,(f) (1) 



o 




Fig. 1. The polyhedron of observing stations (black) and the 
polyhedron of observed bodies (grey). The relative orientation 
of two polyhedrons is estimated from observations of projec- 
tions of vectors between observing stations and observed bodies 
and interpreted as the Earth's rotation. 



where M is the rotation matrix, T'(f) is the translational motion 
of the network of stations, and d{t) is a displacement vector. 
Equations of photon propagation tie the instantaneous vector of 
site coordinates r,(f) with vectors of observed physical bodies 
and their time derivatives. These relationships allow us to build 
a system of equations of conditions. Station position vectors at a 
given epoch and the quantities on the right-hand side of expres- 
sion[T]are estimated from a single least square solution. 

The displacement vector rf,(f) characterizes the motion of an 
individual station, while matrix M and vector T describe the mo- 
tion of the entire network. Assuming the stations are solidly con- 
nected to the Earth's crust, we consider that this part of motion 
represents the motion of the entire Earth. In particular, matrix 
Mit) describes the Earth's rotation. Schematically, the mechan- 
ical model of observations can be viewed as the motion of the 
polyhedron of observing stations with respect to the polyhedron 
of observed bodies (Fig.lTJ. It should be noted that the EOP are 
defined here as the parameters of an estimation model, while in 
the framework of the traditional approach, they are defined as 
angles between big circles on a sphere. 

Since both M(t) and vector d{t) are functions of time, i.e. 
infinite sets of points, they can be evaluated from a finite set 
of observations only in the form of an expansion in some func- 
tions. When we say that the matrix M(t) is determined from ob- 
servations, this should not be understood literally, but instead it 
should be construed that a mathematical model for the depen- 
dence of A1(0 on time is assumed, either explicitly or implicitly. 
The model depends on a finite set of unknown parameters that 
are determined from observations. 

The choice of the mathematical model is not unique. On one 
hand, the mathematical model should approximate the rotation 
with errors comparable to uncertainties of observations during 
the full interval of observations. On the other hand, we should 
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be able to estimate robustly all the parameters of the model. Let 
us consider several approaches. 

2.2. The time series approacii 

The easiest way to represent a rotation is to estimate the matrix 
Ai(t) at certains moment of time and, thus, generate the time se- 
ries. The 3x3 matrix M has 9 elements, but only three of them 
are linearly independent. An arbitrary rotation matrix can be de- 
composed in a product of several elementary rotation matrices 
with respect to coordinate axes at certain angles. Therefore, it is 
sufficient to determine these rotation angles in order to determine 
the matrix A[(t) from observations. 

The fundamental problem is that no observation technique, 
except the laser gyroscope, is sensitive to the instantaneous 
Earth's rotation vector or to its time derivatives directly. The ro- 
tation angles can be derived using the least square estimation 
procedure, together with evaluating other parameters. It requires 
accumulating sufficient amount of data in order to separate vari- 
ables. The estimates of the Earth's rotation angles cannot be 
sampled too fast. A typical sampling rate of estimates is one day, 
since this allows compensation for a certain type of systematic 
error. In some cases the sampling rate can be reduced to several 
hours. 

Unfortunately, one cannot neglect changes in the Earth's ro- 
tation angles during the sampling interval. The accuracy in de- 
termining rotation angles for the 24-hour period is nowadays at 
the level of 2-5 ■ 10 rad. The amplitude of the quasi-diurnal 
motion around axes 1 and 2 is growing with a rate that is an order 
of magnitude of 7 • 10"'^ rad s"'. Therefore, this motion should 
first be separated from the slowly varying components. In the era 
of optical astrometry, some components of this motion, namely 
precession and nutation, were determined separately from ob- 
servations of slowly varying components using a different tech- 
nique and even different instruments. The observations of slowly 
varying constituents in the Earth's rotation angles were corrected 
with a model of the quasi-diurnal motion. iHerring et al.l (Il986l) 
have demonstrated that corrections to the model of the quasi- 
diurnal motion around axes 1 and 2 can be estimated together 
with slowly varying components of the Earth rotation, if the ro- 
tation angles around coordinate axes A,(f) are parameterized in 
the form 

Ai(t) = bi(t) + bi(t)(t - to) + c(t) cos-Q„r + s(t) sin-Q„f 
Mit) = b2(t) + b2(t)(t - to) + c(t) sin-£2„f - s(t) cos-il„f (2) 
A3(f) = b3(t) + bi(t)(t - to) 

where Q„ is the nominal diurnal Earth's rotation rate, 
7.292 115 146 706 707- 10"^ rad s"'. Parameters c{t\ s{t), bi{t) 
are slowly varying functions of time. This approach quickly be- 
came traditional for processing VLBI experiments, and eight 
parameters, bi,b2, bj, bi,b2, bj, c, s are routinely determined for 
each individual 24 hour observing session. 

2.3. Limitations of tiie time series approacii 

However, it is important to realize the limitations of the time se- 
ries approach. First, the raw time series of estimates provides the 
values of rotation angles only at specific discrete moments. They 
do not determine a functional dependence of rotation angles on 
time. An end user needs to have a tool for computing Earth's 
orientation at any moment of time within the interval of obser- 
vations. Thus, the raw time series are the basis for the second 
step of processing that involves smoothing and interpolation. 



Smoothing and interpolation of the time series q, s^, bik, b2k, b^k 
implicitly assumes that A,(0 satisfies some mathematical model 
that appears to be different from the model in expression|2]used 
in the estimation process. The resulting smooth function of rota- 
tion angles does not provide the best fit to observations; if it did, 
no smoothing would have been needed. 

Second, at the present level for accuracy of observations, 
the estimation model corresponding to Eqs. |2]is not adequate: 
one cannot neglect changes in c(t), s(t) and bi(t) over the inter- 
val of estimation, typically 24 hours. Adjusting time derivatives 
c(t), i(f), bi(t) makes estimates of these parameters so strongly 
correlated that they do not have a practical value. 

Changes in the a priori model for slowly varying compo- 
nents of rotation angles affect all estimated parameters, includ- 
ing c(t) and s{t). In order to demonstrate this, two VLBI so- 
lutions using 3563 twenty four hour observing sessions from 
1984 through 2006 were computed. The USNO Finals EOP 
time series of pole coordinates and UTl-TAI with a time span 
of 1 day (Dick and Richter 20040 were used as the a priori 
model in the reference solution. The Gaussian noise with stan- 
dard deviation 1 nrad was added to all components of the USNO 
Finals EOP series, and these modified time series were used 
in the trial solution. The rms of differences in the total val- 
ues of b{t), c{t), s{t), i.e., the sum of a priori values and adjust- 
ments over the 24 hour time intervals, was 0.14 nrad for b(t) 
and 0.16 nrad for c(t)ands(t). Since the accuracy of estimates of 
b(t), c(t), s(t) from 24 hour VLBI experiments is cuiTently at a 
0.3 nrad level, the accuracy of the a priori EOP series should be 
better than 1 nrad in order to reduce the contribution its errors to 
a negligible level: 1/2 of the random error in estimates. In a sim- 
ilar way, the change in the a priori model for the quasi-diurnal 
motion also affects estimates of c{t), s(t) and bi{t). Although one 
can expect that a continuous process or refining the a priori 
model and subsequent least square estimation should converge, 
this does not happen in practice. It is known among analysts who 
process raw data that, if the initial a priori values are changed, 
the total angles, i.e. the sum of the a priori and the adjustments, 
come out different. Researchers who process time series are not 
always aware of these complications and tend to consider the re- 
sults of processing the same observations by different analysts 
as independent "data", so they attribute the differences between 
them to so-called "analyst noise". These discrepancies occur due 
to an internal inconsistency between the estimation model, the 
a priori model, and the post-processing procedure. 

Third, the second step in the analysis, smoothing and inter- 
polation, is rather subjective. A different degree of smoothing 
produces a different series. 

Finally, the time series cannot be used directly for making 
an inference about the physical processes that affect the Earth's 
rotation. The time series are transformed by various analysis pro- 
cedures. The dependence of the series on the a priori model and 
the correlations between the elements of time series are usually 
ignored. The correlations between the elements are not zero, 
because the elements themselves were estimated together with 
other parameters like global site velocity or source positions that 
are considered common for the entire interval of observations. 
Although these coiTelations are not strong, typically at a level 
of 0.1, their contribution is significant when long time series are 
processed. 

Due to the complexity of the a priori model, analysts who 
process the time series of estimates usually do not handle the to- 
tal angles of the Earth's rotation, but rather adjustments to the 

^ Available on the Internet at ftp : //maia.usno.navy.mil/ser7/finals . ^ 
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a priori values, tacitly assuming that analysts who processed the 
raw data strictly followed a standardized procedure for data re- 
duction. In reality this is often not the case. This creates an ad- 
ditional source of confusion and errors. 

These complications prompt us to look for a one-step proce- 
dure of estimation of the Earth's orientation parameters. 



3. Representation of the Earth rotation in the form 
of the expansion into basis functions 

While Eqs. |2]represent rotation angles within a short period of 
time, 24 hours, they are not adequate for a longer period of time. 
We need to find a mathematical model which would be valid for 
the entire period of observations, i.e. several decades. The ma- 
trix Ma has a non-linear dependence on its arguments. A linear 
estimator, the least square method, allows us to evaluate not the 
matrix itself, but its small perturbation. The coordinate transfor- 
mation of a vector r from the terrestrial coordinate system to the 
celestial coordinate system is then written as 



rc = Ma(t) + [qeit) + qAt)) x 



(3) 



where and designate the coordinates of the vector r in the 
celestial and terrestrial coordinate systems, respectively, q^it) is 
the vector of a small perturbational rotation, qa{t) is the small 
a priori rotation vector in the terrestrial coordinate system, and 
Mait) is the a priori matrix of finite rotation. Vectors qe{t) and 
qaif) are small in the sense that we can neglect squares of their 
components. The vector qa{t) can be set to zero by an appropri- 
ate choice of the matrix Maif). Considering that the accuracy of 
determination of rotation angles averaged over a 24 hour period 
is currently at the level of 3 -lO"'" rad, and the accuracy of esti- 
mates of amplitudes of harmonic constituents averaged over the 
period of 20 years is at the level of 10 " rad, the components of 
vectors qa{t),qe{t) should not exceed 3 ■ 10"^ rad. It should be 
noted that these requirements on accuracy of the a priori model 
are three orders of magnitude weaker than those needed for an 
unbiased estimation of time series. 

In order to find an appropriate basis for expanding of qAt), 
we need to use an a priori knowledge of the process under con- 
sideration. The Earth's rotation can be considered in terms of a 
response to external forces. The external forces that affect ro- 
tation of the solid Earth are caused 1) by redistribution of geo- 
physical fluids and 2) by tidal attraction of external bodies. The 
first process is not predictable and is dominating at frequencies 
by modulo much less than the diurnal frequency Q„. The tide- 
generating potential exerted by external bodies can be consid- 
ered to be known precisely. Its spectrum has a comb of very 
sharp lines as shown in Fig.|2l 

To characterize the Earth's response, we should take into ac- 
count that the triaxiality of the Earth's inertia tensor (B - A)/C 
is small, about 2 ■ 10"^. Therefore, the differential equations of 
the Earth's rotation are linear First, this leads to decoupling ro- 
tation around the axes 1 and 2, i.e. the polar motion, and rotation 
around the axis 3, the diurnal motion. Second, the response to 
harmonic external forces will result in harmonic variations of 
the component 2 of q,, with the same amplitude as component 
1 with the phase shifted by -;7r/2. Third, the excitation at the 
diurnal frequency will result in the appearance of cross-terms 
f sin-Q„ and f cos-Q„ ( lMoritzllI987h . 
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Fig. 2. The logarithm of the power spectrum of the tide generat- 
ing potential in kg^m'^s as a function of the angular frequency. 



Considering these factors, the following mathematical model 
for the the vector of a small perturbational rotation is proposed: 

/ n-l N 

Yj fikBl'(t) +^(P5cosw,„f + P'jsinajjt) 



qeit) 



k=l-m 



+ t (S'' cos-f2„ t + S ' sin-Q„ f) 



2 f2kB'^(t) +Yj{P'jSmcjjt - P'j cos cj J t) 



+ t [S " sin-Q„ t - 5*cos-Q„f) 



(4) 



f,k B"^{t) + ^ [e) cos ojj t + E] sin ojj t) 

\ k=l-m J=l 

where B™(0 is the B-spline function of degree m deter- 
mined at a sequence of knots fi_„,, f2-,„, ■ • ■ , fo, fi, ■ ■ ■ f*; 
ojj are the frequencies of external forces; the coefficients 
fill, P'j, P'j, S'^,S'',E'j,E'j; are the parameters of the expansion, 
Q.„ is the nominal frequency of the Earth's rotation. Here n is 
the dimension of the B-spline basis and is the dimension of 
the Fourier basis. Thus, the vector of small perturbational rota- 
tion is expanded into the basis of B-splines, which is orthogonal 
over the entire period of observations, and the basis of harmonic 
functions, which is orthogonal in the range (-00,-1-00). The first 
basis approximates the slowly varying component in the Earth's 
rotation, the second basis approximates the quasi-diurnal com- 
ponent, as well as other harmonic constituents of the Earth's ro- 
tation. 



3.1. The B-spline basis 

The B -spline basis functions were first introduced by lSchonbergI 



(I1946). The B-spline function of degree m depends on time and 
on a monotonically nondecreasing sequence of n knots at the 
interval [ti,t„]. In order to introduce splines, let us extend this 
sequence by adding m elements at the beginning of the sequence 
and m-l elements at the end of the sequence such that fi_„, = 
ti-m = . . . - to - ti and t„ = f„+i - t„+2 - . .. - tn+m-\- At a 
given extended sequence of n H- 2m - 1 knots, n + m-l B-spline 
functions with the pivot element k e \ - m, 2 - m, ...n-l are 
defined through a recursive relationship. 

The B-spline of the 0-th degree with the pivot knot k e [ 1 , n- 
1] on the knots sequence (fi, f2, ... ,f„), such that ti < t2 < 
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Fig. 3. The logarithm of the power spectrum of quasi-diurnal 
variations in ^i, ^2 components of the rotation vector according 
to the REN-2000 expansion in rad^ as a function of the angular 
frequency. 



< tk, is determined by 

1, if f e {tk,h+i) 



0, otherwise 



(5) 



The B-spline of the mth degree with the pivot knot k e [ \-m, rt-l ] 
on the extended sequence of knots (fi_,„, f2-m, ■■■ tn+m-i, fn+m-i) 
is expressed via the B-splines of the m- 1th degree as 



t - tk 

Bm - : — 



^k+m ^k 



(6) 



Computation of B-splines is as simple as computation of 
other polynomials. Similar simple recursive relationships exist 
for derivatives of B-splines and integrals. The B-spline of de- 
gree m with the pivot element k is non-zero only at the interval 
{tk, tk+m+\)- It can be proved that a sequence of n + m - 1 B- 
spUne functions of degree m with pivot elements k e I - m, 2 - 
m, . . . ,n - 1 forms a basis on the interval [fi, f„]. The proof of 
this and many other useful theorems related to B-splines can be 
found in Niirnberger (1989.). 

In general, knots can be selected arbitrarily. Test runs have 
shown that a set of B-spline functions of the 3rd degree with 
equidistant knots with a time span of 3 days for components 1,2, 
and 1 day for component 3 of the vector ^p(f) adequately repre- 
sents the slowly varying component of the Earth's rotation. Weak 
constraints on values of B-splines, its first and second derivatives 
can be imposed to ensure the stability of the solution at intervals 
with considerable gaps in observations and at the beginning and 
the end of the data set. 

3.2. The Fourier basis 

JVIodeling the quasi-diurnal components is more challenging. 
The tides exerted by the Moon and the Sun cause variations in 
sea currents and the sea surface at tidal frequencies. These vari- 
ations excite changes in all components of the Earth's rotation. 
The resonance near the retrograde diurnal frequency causes a 
significant amplification at that frequency band for components 
1 and 2 of the vector qdl). The theoretica l spectrum of this mo- 
tion referred to as nutation computed by ISouchav & Kinoshital 
([1996, 1997) and Souchay et al. ( 1999^) for the model of the rigid 
Earth, the REN-2000 expansion, is presented in Fig. [3] 

The problem is that the spectrum is very dense, and observa- 
tions during a finite period of time cannot resolve all the con- 
stituents. The REN-2000 spectrum has 560 constituents with 



amplitudes greater than 10"" rad and 1551 constituents with 
amplitudes greater than 10 rad with the frequency difference 
between some of them as low as 10"'^ rad s"'. 

Several strategies can be used for overcoming this problem. 
First, we can select frequencies with maximal amplitudes from 
the theoretical spectrum and ignore constituents with an angular 
frequency separation less than w„„„ = Itt/AT, where AT is the 
interval of observations. No signal will be mismodeled if the 
frequency separation between the constituents Ao) « a>,nw, since 
in this case the two constituents will be indistinguishable. 

However, if the constituents are not very close, the mismod- 
eled signal will leak into adjustments at other frequencies. The 
sidelobe with the amplitude A2 and frequency a>2 of the main 
constituent with the frequency a>i can be omitted if the quan- 

tity A2—7 is less than a certain threshold. Depending on 

the threshold level, there are several hundred constituents in the 
tidal spectrum for which this condition is not valid. 

One way to mitigate this problem is to estimate the ampli- 
tude of close sidelobes, together with the amplitude of main 
constituents, and to impose strong constraints on the amplitude 
of sidelobes by using some a priori information. It is plausible 
to assume that the a posteriori amplitudes of constituents of the 
quasi-diurnal motion for the real Earth differ from the theoretical 
amplitudes computed for the rigid Earth by multiplicative factors 
called transfer function, which is a smooth function of frequency 
according to theory. Then we can assume that the transfer func- 
tion for two close constituents with theoretical amplitudes P be 
the same, i.e. the ratio of complex amplitudes A of two close 
constituents is the same as for the a priori rigid Earth amplitudes, 
and, therefore, should satisfy this equation: 



P'l+iPl 



•/Af 



iPl 



■iA' 



(7) 



Although this approach reduces the leakage from a mismod- 
eled signal, it is not fully satisfactory. In general, using strong 
constraints is undesirable, since this introduces a bias in esti- 
mates. The validity of Eq.|7]cannot be confirmed or refuted from 
observations. It comes from a theory. But if the estimation model 
implicitly incorporates theoretical assumptions, strictly speak- 
ing the estimates cannot be used for validation of the theory. 
Although Eq. Q for EOF variations caused by the tidal poten- 
tial exerted by external bodies has a sound theoretical basis, we 
should bear in mind that the ultimate goal of comparing theo- 
retical predictions with observations is to check the validity of 
assumptions built into the foundation of the theory and to make 
a judgment whether the model is complete or not. If there are un- 
accounted additive constituents at these frequencies, for exam- 
ple, caused by the free motion, by the atmospheric, or by oceanic 
excitation, the Eq. |7]may not be valid. 

An alternative to constraining sidelobes is the strategy of es- 
timating a wide range of constituents that are multiples of w,„,„ 
or, in other words, indirectly performing the discrete Fourier 
transform of the perturbational rotation. With this approach, in 
general we are in a position to discard our a priori knowledge 
about the frequency structure of the signal. Estimating the signal 
at discrete frequencies that are multiples of comin, from the zero 
frequency through the Nyquist frequency, recovers any signal ac- 
cording to the sampling theorem. However, this kind of approach 
applied to estimating the vector qe(t) has a practical value only 
if the number of non-negligible constituents in the discrete spec- 
trum is significantly less than the total number of samples. Since 
the spectrum of the tide-generating potential consists of a set 
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of discrete frequencies that are not commensurate to each other, 
the frequencies that are multiples of (y„„„ cannot coincide with 
all tidal frequencies. If the amplitude of a narrow-band harmonic 
signal is not estimated at its frequency, but estimated at a set of 
nearby frequencies that are multiples of w„„„, the signal will be 
recovered only partially. The wider the range of frequencies, the 
better the approximation. The rate of convergence depends on 
the amplitude of the signal and the difference between its fre- 
quency and the closest frequency used for estimation. Selection 
of reasonably good a priori qa{t) values may significantly reduce 
the number of frequencies needed for estimation to reach a given 
level of accuracy of approximation. 

Other important constituents of the signa l at the retrog rade 
diurnal band are the free ne ar-diurnal wobble (Moritz 1987) and 
the at mospheric nutations CBizouard et al. 1998; Yseboo dt et aTl 
I2OO2I) . Since this signal is excited by a broad-band stochastic 
process, it is expected that these constituents in the Earth's ro- 
tation are also relatively broad-band. To model this signal, the 
constituents at frequencies within the range of that band need 
to be added to the list of constituents at tidal frequenci es. It 
follows from the sampling theorem of iKotelnikovl (Il933l) that 
a band limited signal with frequencies in the range of [w/, w/,] 
is completely recovered when the estimates of the sine and co- 
sine amplitudes of the spectrum are made at discrete frequencies 

[(Dl, (L>l + (L>„,i„, 0)1+ 2 (L>„i„, ... (Dl + (N -I) (Jj„u„, (jJh] . 

The tidal spectrum also has constituents with low frequen- 
cies, so-called zonal tides. They affect component 3 in the vector 
of the perturbational rotation. Their contribution dominates the 
rate of change for this component. It would be desirable to esti- 
mate the complex amplitude of this variations. Since the residual 
rate of change of is a factor of 3-10 less, constraints on a rate 
of change for the residual component of qe, modeled with an 
expansion over the B-spline basis can be set stronger without 
introducing a bias in the estimates. This improves the solution 
stability during intervals of time with gaps in observations. For 
the same reason, it would be desirable to estimate variations in 
components 1 and 2 of the Earth's rotation vector at the annual 
and Chandler frequencies: 1.990968 ■ 10"^ and 1.678 • 10"^ rad 
s"' , respectively. 



3.3. Decorrelation constraints 

It should be noted that the estimates of harmonic constituents 
with lower periods than the time span between nodes of B-spline 
will so highly correlate with B-spline coefficients that the sys- 
tem of equations will be very close to singular. Decorrelation 
constraints on coefficients of the B-spline should be imposed in 
order to overcome this problem. We require that the product of 
expansion over basic B-spline and Fourier functions for the jth 
frequency be zero over the interval of observations: 



to ^*=1"'" 

'r( "I 



\k= 1 -m 



dt^O 



dt^O 



(8) 



This is reduced to 

, +00 



fk [ Bl^it) cos LOj t dt ^0 

=1-'" -"i 

, +00 

n-l p 

^ A I B™(0sinw; ft/f = 



(9) 



k=\-m 



Thus, two constraint equations for each frequency are to be 
imposed. The Fourier integral of a B-spline of the mth degree in 
Eq.|9]on a knots sequence (t^, ti^+i, ... ,t„) with the pivot knot k 
such that A; - m < n - 1 is expressed through the Fourier integral 
of a B-spline of the m - 1 th degree: 



I 



to 



«^(0e dt = \Bi^(t„)e "-BAtije ' + 



(10) 



to{tk+m - tk) * 



<y(ft+I - tk+m+l) 



The Fourier integral of a B-spline of the 0-th degree on the 
same sequence (tk, tk+\, ... , f„) with the pivot knot / is 



+00 



I B'l{t)e^'''dt^ - 
(jj 

—00 



(11) 



4. Analysis of VLBI observations 

4.1. The VLBIdataset 

A set of estimates of group delays at frequency bands cen- 
tered around 2.2 and 8.6 GHz from January 1984 through 
January 2006 was used to validate the proposed approaches. The 
International VLBI Service for Geodesy and Astrometry (IVS) 
(Vandenberg 1999) provides online access to the collection of 
all observations made in the geodetic mode under various as- 
trometric and geodynamics programs from 1979 through now at 
httpT77ivscc.gsfc.nasa.gov The VLBI data set shows a 
substantial spatial and time inhomogeneity. Typically, observa- 
tions are made in sessions with a duration of about 24 hours. 
Observations were sporadic in the early 80s, but in January 1984 
a regular VLBI campaign for the determination of the Earth ori- 
entation parameters started first with 5-day intervals, from May 
1993 with weekly intervals, and from 1997 twice a week. In ad- 
dition to these observations, various other observing campaigns 
were running. On average, 150 sessions per year have been ob- 
served since 1984. 

During that period 153 stations participated in observations, 
although a majority of them observed only during short cam- 
paigns. The observations at stations that participated in less than 
20 000 observations, and the stations that only participated in 
at regional networks with sizes of 2000 km and less were dis- 
carded. Forty four stations remained. Observations of sources 
that were observed in less than 4 sessions and gave less than 64 
usable pairs of dual-band group delays were excluded. The data 
before January 1984 were also discarded. In total, ~5% of the 
observations were excluded, and the remaining data from 3563 
sessions between January 1984 to August 2006, more than 4.6 
million of dual-band pairs of group delays, were used in the anal- 
ysis. 

The number of participating stations in each individual ses- 
sion varies from 2 to 20, although 4-7 is a typical number. No 
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station participated in all sessions, but every station participated 
in sessions with many different networks. All networks have 
common nodes and, are therefore, tied together. Networks vary 
significantly, but more than 70% of them have a size exceeding 
the Earth's radius. 



4.2. Theoretical model 



The state of the art theoretical models were used for computing 
the theoretical time delay and its partial derivatives . The proce- 
dure in general follows the approach presented by Sovers et al.l 
( 119981) with some minor refinements. The expr ession for time 
delay derived by iKopeikin and Schafeii (Il999l) in the frame- 
work of general relativity was used. The displacement caused 
by the Earth' s tides were computed using a rigorous algorithm 
iPetrov & Mai (|2003|) with a truncation at a level of 0.05 mm 
using the nu merical va l ues of the generalized Love numbers 
presented by iMathewsl (l200lh . The displacements caused by 
ocean loading were computed by convolving the Greens' func- 
tions with ocean tide models using the NLOADF algorithm of 
lAgnewl (Il997ii) . The GOTOO model ( Ravlll999l) of diurnal and 
semi- diurnal ocean tides, the NA099 model jMatsumoto et al.l 
I2OOOI) of ocean zonal tides, the equilibrium model of the pole 
tide and the tide with period of 18.6 years were used. The at- 
mospheric pressure loading was computed by convolving the 
Greens' functions with the numerical model of the atmosphere 
NCEP Reanalysis jKalnav et al.lll996(). The algorithm of c om- 



used: 



putations is described in full details in 



Petrov & Bovl (120041) . 



The a priori path delay in the atmosphere caused by the hy- 
drostatic component was calculated as a product of the zenith 
path delay computed on the basis of surface pressure using 
the Saastamoinen (1972 ) expr ession and the isobaric mapping 
function (iNiell & Petrovll2004t) computed using the geopotential 
height of the 20 kPa pressure layer provided by the numerical 
weather model NCEP Reanalysis. The isobaric mapping func- 
tion describes the dependence of path delay on the angle between 
the local axis of symmetry of the atmosphere and the direction to 
the observed sources. The direction of this axis from the zenith 
was considered to coincide with the normal to the surface of the 
geopotential height at the 20 kPa pressure level. This normal was 
computed using the NCEP Reanalysis dataset. 

Since the accuracy requirements to the of the a priori Earth 
rotation model are very low in the framework of the present ap- 
proach, we can exploit this to use the simplest possible model. 
The following expression for the a priori matrix of the Earth's 
rotation M„{t) according to the Newcomb-Andoyer formalism 
was used: 



Ma(f) = ^3(^0) ■ Kii-Oa) ■ K^iz) ■ Kii-eo) ■ H^i^^) ■ 
??i(eo + Ae)-??3(-5) 



(12) 



9o = 9oQ + Omt + 9o2 f 

Z = Zo + Zl t + Z2 



(13) 



where Hi is a rotation matrix around the axis /. For the variables 
^0, ^0, z, eo, Aiff, Aeo, S, the following simplified expressions were 



fo - foo + £01 f + £02 1 

2 

AiA = ^Pj sin (aj+/3jt) 
i 

Ae = ^^Cj cos (a j + fljt) 

S = So+7t-Eo + (n„+(oi+Zl-E{}t + {(02+Z2-E2)f 
2 

+ At// COS 60 - ^ {E^ COS y, t + E- sin 7, t) 



Here f is TAI time elapsed from 2000 January 01, 12 hours. It 
should be noted that expi'essions[T3]diffe rs from those used in the 
framework of the traditional approach dJMcCarthv & Petit. eds.l 
2004). Some of these parameters were taken from theory, some 
of them were found with the LSQ fit of time series of adjust- 
ments of pole coordinates and UTl-TAI. The numerical values 
of parameters used in data analysis are presented in the online 
Tab. [3] The rms of the adjustments of the perturbational vector 
of the Earth rotation with respect to the a priori matrix presented 
in expressions [T2l and [T3] over the period of 1984-2006 was less 
than 2.0 ■ 10"* rad for each component. 

4.3. Basic estimation model 

Several solutions were produced. Each solution used the basic 
parameterization which was common for all runs, and a specific 
parameterization for an individual solution. Basic parameters be- 
long to one of the three groups: 

— global (over the entire data set): positions of 598 sources and 
proper motions of 169 sources; positions and velocities of 
44 stations; coefficients for the expansion into B-spline ba- 
sis of positions of 7 stations, DSS15, DSS65, GILC REEK 
HRAS.085, MEDICINA, MOJAVE12, PIETOWN, ( Petrovl 
200_5); coefficients for harmonic position variations of 
all sites at the an nual and the semi-annual frequency 
jPetrov & Mall2003h ; coefficients of the B-spline for mod- 
eUng the perturbational Earth's rotation vector qeit) at a set 
of equidistant knots with a time span of 3 days for compo- 
nents 1 and 2 and with a time span of 1 day for component 
3. 

— local (over each session): tilts of the local symmetric axis 
of the atmosphere for all stations and their rates, station- 
dependent clock functions modeled by second order poly- 
nomials, baseline-dependent clock offsets. 

— iegwenfeii (over 0.33-1.0 hours): coefficients of linear spline 
that models atmospheric path delay (0.33 hour segment) and 
clock function (1 hour segment) for each station. The es- 
timates of clock function absorb uncalibrated instrumental 
delays in the data acquisition system. 

The rate of change for the atmospheric path delay and clock 
function between two adjacent segments was constrained to zero 
with weights reciprocal to 1.1 ■ 10"'"* and 2- lO"''*, respec- 
tively, in order to stabilize solutions. Strong no-net-translation 
and no-net-rotation constraints were imposed on the adjustments 
of site positions and velocities, as well as no-net-rotation con- 
straints were imposed on adjustments of source positions and 
proper motions, in order to solve the LSQ problem of incom- 
plete rank. Weak stabilizing constraints were imposed on B- 
spline coefficients to constrain to zero qe{t), its first and second 
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derivatives at each knot. The reciprocal weights of constraints 
for the values of qeit), components 1 and 2 were 5 ■ 10"^ rad; 
the reciprocal weights for the first time derivatives were 
5 • 10"'"* rad s and for the 2nd derivatives 3 ■ 10 '"^ rad s"^. 
The reciprocal weights in qe{t), qdt), qe{t) for component 3 were 
5 • 10-^ rad, 3 ■ lO^'^ rad s"', and 6 ■ IQ-^'' rad s-^. 

4.4. VLBI solutions 

Several solutions have been computed with different a priori vec- 
tors of the perturbational rotation qa(t) and with a different set 
of harmonic constituents of the vector qdt). 

In solution A, the a priori vector q„{t) was set to zero. A 
set of sine and cosine amplitudes of the perturbational rota- 
tion vector was estimated as global parameters. The frequen- 
cies of these estimates were selected according to the follow- 
ing process. First, the frequencies for components 1 and 2 of 
this vector in the range of the near-diurnal retrograde wobble 
[-7.310955 ■ 10-^ -7.298755 ■ 10"^ rads"'] a nd hypothetica l 
near-diurnal prograde wobble (see, for instance. lOehantl i 1 993b ) 
[-7.284405 ■ 10-^ -7.273275 ■ 10"^ rads"'] were sampled with 
step 9.0 ■ 10"^ rads '. This step of frequency sequences corre- 
sponds to 2n/AT where AT" is the interval of observations, 22.6 
years. Then the frequencies of constituents with the amplitudes 
exceeding 1 ■ 10 " rad in the REN-2000 expansion were added 
to the list of constituents in order of decreasing their ampli- 
tudes, provided the minimal difference in the frequencies of the 
constituents was less than 9.0 ■ 10"^ rads If two constituents 
have difference in frequencies less than that, the constituent with 
smaller amplitude was discarded. 

Second, all constituents at positive and negative diurnal, 
semi-diurnal, and ter-diurnal bands of the tide-generating poten- 
tial with amplitudes greater than 0.003 of the amplitude of the 
M2 tide were selected. At negative frequencies, all three compo- 
nents of the vector q^it) were estimated, while only the compo- 
nents 1 and 2 were estimated at positive frequencies. The same 
rejection criteria for the constituents with close frequencies was 
enforced. 

Third, the harmonic signal in components 1 and 2 of q^ at the 
Chandler and annual frequencies, and the harmonic signal in the 
component 3 at 14 frequencies of zonal tides with amplitudes 
greater than 3% of the amplitude of the tide generating poten- 
tial at the frequency 5.323414 ■ lO"*" rad s"' were estimated. 
Decorrelation constraints were imposed on estimates of B-spline 
coefficients. 

Solution B is similar to solution A, but constituents at the 
76 frequencies that have a close companion, a sidelobe, and that 

satisfy the criteria in Sect. [3] — — > 1 ■ 10 " rad, were 

not rejected as in solution A, but remained on the list. Strong 
constraints with reciprocal weig ht 10 ^^rad^ in the form of Eq.Q 
were imposed. The purpose of this solution was to investigate the 
effects of omitted sidelobes. 

A family of solutions C was computed. The purpose of this 
solution was to evaluate a harmonic signal at non-tidal frequen- 
cies. In addition to the frequencies used in the solution A, a pos- 
sible non-tidal signal was sought in eight frequency bands at the 
frequencies equally sampled with the step 9.0 ■ 10"^ rad s ' . The 
low and high edges of each frequency band are shown in Tab.[T] 
The total number of constituents of the perturbational rotation 
vector estimated for components 1 and 2 was 20 093 and and for 
component 3 was 9885. Since the total number of global param- 
eters was at a level of 70 000, well beyond the capabilities of a 



Table 1. The range of frequency bands for estimation of a non- 
tidal signal. The last column refers to components of q^it) vector 
of perturbational rotation that were estimated. 



low 






Frequency band 
high 








V-UllipUllCllLS 


-2.95 • 


10- 


-4 


rads"' 


-2.85 ■ 


10- 


-4 


rad s"' 


1,2,3 


-2.27 • 


10- 


-4 


rad s"' 


-2.07 ■ 


10- 


-4 


rad s"' 


1,2,3 


-1.60- 


10 


-4 


rads"' 


-1.32- 


10 


-4 


rad s"' 


1,2,3 


-0.97 • 


10- 


-4 


rad s"' 


-0.52 • 


10 


-4 


rad s"' 


1,2,3 


0.52' 


■ 10 


-4 


rad s"' 


0.97- 


10 


-4 


rad s"' 


1,2 


1.32' 


■ 10 


-4 


rad s"' 


1.60- 


10 


-4 


rad s"' 


1,2 


2.07 ■ 


■ 10 


-4 


rad s"' 


2.27 ■ 


10 


-4 


rad s"' 


1,2 


2.85 ■ 


10 


-4 


rad s"' 


2.95- 


10 


-4 


rad s"' 


1,2 



personal computer, 20 individual solutions were performed. In 
each individual solution of family C, the sine and cosine ampli- 
tudes of constituents at ~ 2000 frequencies were estimated. The 
estimates of constituents are correlated at a level of 0.02-0.20, 
since other common parameters were estimated, such as source 
coordinates and station positions. Strictly speaking, the proce- 
dure of estimating the spectrum by parts is not perfectly correct. 
But it was assumed that such a procedure may result in a false 
detection, but not miss the signal present in the data. In the fi- 
nal run of a solution of family C, 680 non-tidal frequencies were 
selected. The signal was detected at a 95% confidence level at 
these frequencies in the previous runs. 

The purpose of solution D was to investigate whether the 
process of frequency selection for solution A picked up all 
the signals in the diurnal band. In this solution the a pri- 
ori vector qa{t) was generated from the REN-2000 nuta- 
tion expansion. It had all the constituents with amplitudes 
greater than 1 ■ 10"^ rad, except the constituents with frequen- 
cies -7.331937- 10'^ -7.293186- 10"^ -7.291047- 10"^ and 
-7.252295- 10"^ rad s"', because they had already been in- 
cluded in the matrix Ma(t). For selection of the frequencies, at 
which the harmonics of the perturbational rotation vector was 
to be estimated, the time series of q^it) + qgif) - q'^it) on the 
interval [1984.0, 2006.6] with a step of 0.125 days were pro- 
duced, and its power spectrum was computed. Here q^(t) is the 
vector of the adjustments of perturbational rotation in solution 
B and vectors q^, where are the a priori perturbational vec- 
tors used. The frequencies of constituents run from through 
2.9 • lO""* rads-i with step 9.0 • 10"'' rads"^ Not all of them 
were estimated, but only those that had the frequency by module 
less than 9.3 ■ 10"^ rad s"' and that satisfied one of the three con- 
ditions: a) to be in the range [7.20 - 10"^ 7.38 - 10"^] rad s"' for 
taking the near-diurnal wobble and the signal excited by the at- 
mosphere into account, or b) to have the square root of the power 
greater than 110"" rad for taking the signal at tidal frequencies 
into account, or or c) to be in the range of non-tidal frequencies 
identified in solutions C. In total, 3676 parameters at 1706 fre- 
quencies, which are a multiple of the frequency 9.0 ■ 10"^rad s"' , 
were selected. 

The software Calc/Solve was used for these solu- 
tions. Results of analysis are available on the Web at 
http://vlbi.gsfc.nasa.gov/ermi 

5. Discussion of results 

5. 1 . Differences with respect to the traditional approach 

Analysis of the results showed that the approach for the direct 
estimation of the perturbational rotation vector in the form of 
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expansion into basic functions from VLBI observations allows 
us to represent the Earth's rotation adequately. The fit of the least 
square solutions, 21.9 ps, was the same as the fit in solutions 
that followed the traditional time series approach. The slowly 
varying component was compared with the USNO EOP model. 
The USNO Finals EOP time series with 1 day steps was derived 
by averaging results of analysis of GPS and VLBI observations 
and smoothing. To compare results, the USNO series of polar 



motion, Xp, Yp, and UTl-TAI with respect to the a priori Earth % 



rotation model in expressions[T3] 

q'i = Yp(t) 
ql = Xp{t) 

q"^ = k{UT\ - TAI){t) + (£0 + £1 f + E2 fi) ■ 



^ (fi^ cos ji t + Ej sin ji tj + J(i}f + Ai}f) Ae sin eo dt 



(14) 



where k - -(Q„ + zoi + ^01) ■ 86400/27r, and parame- 
ters -Eo, -El , E2, E":, Ep Ji, i//, Ai//, Ae, eo, £2„, z, ( are from expres- 
sions [13] The coefficients of the interpolation spline for q" were 
computed. These coefficients form the USNO Earth rotation 
model as a continuous function of time. Since the GPS results 
almost entirely dominate components 1 and 2 of the Earth rota- 
tion vector from that model, they can be considered independent 
from our analysis of VLBI observations. 

The differences for component 1 between the USNO model 
and our results from solution B after removal the contribution 
of harmonic variations with periods less than 2 days are shown 
in Fig. |4] No pattern of systematic differences is revealed. The 
statistics of these differences for all three components of the 
small vector of the Earth rotation and their time derivatives com- 
puted at the equidistant grid with time interval 2.5 hours are pre- 
sented in the 1 st row of Tab. |2] 

Since the VLBI observations are not carried out continuously 
due to budget limitations, the accuracy of the Earth orientation 
model is the highest within an interval of observations and the 
lowest at moments of time when there were no observations. In 
the framework of the traditional approach, the EOP are estimated 
on moments of time in the middle of a 24 hour observing session. 
The statistics of the differences of the EOP series from analysis 
of VLBI observations gs f2006c|] for moments in the middle of 
1426 observing sessions are shown in the 3rd row of Tab.|2l For 
comparison, the EOP were computed from results of solution B 
at exactly the same epochs, and these statistics of the differences 
with respect to the USNO model are presented in the 2nd row of 
this table. 

Analysis of statistics shows that the differences in compo- 
nents 1 and 2 of the Earth's orientation according to the pro- 
posed and traditional approaches do not exceed 20%. At the 
same time the proposed approach gives the estimates of all the 
components of the Earth's angular velocity vector by a factor of 
1.5-2.0 closer to the GPS results than the estimates produced in 
the framework of the traditional approach. According to the tra- 
ditional approach, the EOP rates and nutation daily offsets are 
computed for each session independently, which makes them 
less stable. With the proposed approach, at a given epoch several 
experiments contribute to estimates of EOP rate, which makes 
them more robust. 

Analysis of the differences in amplitudes of the harmonic 
terms of components 1 and 2 of the vector of perturbational ro- 
tation at the retrograde diurnal band with respect to the semi- 
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Fig. 4. The component 1 of the residual perturbational vector 
with respect to the Earth rotation USNO Finals EOP model. 
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Fig. 5. The time series of the estimates of the daily offsets of 
nutation in obliquity when the empirical Earth rotation model 
from solution B was used as the a priori. The wrms is 3.9 ■ 10 
rad. 



empirical MHB2000 expansion (iMathews et al.ll2002h . showed 
they can reach 0.2 nrad for some terms. Detailed analysis of 
these differences is beyond the scope of the present paper In 
order to test results, the empirical Earth rotation model from so- 
lution B was used as the a priori for the solution that estimated 
the time series of daily offset to nutations. The weighted root 
mean square of the differences for the period of [1996.0, 2006.0] 
is 0.39 nrad when results of solution B were used as the a priori, 
and 0.98 nrad when the MHB2000 was used. The daily offsets to 
nutation in obliquity Ae(f) with respect to both models are shown 
in Figs. 



Available on the Web at |http://vlbi.gsfc.nasa.gov/solutions/2006c| 



5.2. Harmonic components in the Earth's rotation 

Analysis of estimates of the harmonic components showed ex- 
cessive power near the frequency of the near-diurnal retrograde 
wobble, as was expected. The spectrum turned out rather broad, 
spanning a rather wide band, and it partly overlaps with the tidal 
frequency -7.312026 ■ 10"^ rad s ' that corresponds to the an- 
nual retrograde nu tation as shown in Fig. [7] It was found by 
iHerring et alj ( 11986 ) that the near-diurnal wobble cannot be rep- 
resented by a purely harmonic model with a constant amplitude. 
This means that when this component of the Earth's rotation is 
represented in the frequency domain, several constituents in the 
spectrum will correspond to it. 

Analysis of the results of the C family solutions revealed sev- 
eral constituents with the non-tidal signal. The spectrum of com- 
ponents 1 and 2 in the vicinity of the frequency -2Q„, i.e. the K2 
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Table 2. The weighted root mean squares of the differences between estimates of the Earth rotation model from analysis of VLBI 
observations and the USNO Finals EOF model for the period of [1996.0, 2006.0]. The statistics in rows 1 and 2 correspond to 
solution B, which follows the proposed approach. The statistics in row 3 correspond to solution gsf2006c, which follows the 
traditional approach. 



Solution q\ q2 qi qi q^ 

ERM all a79~l(P7ad a99~l(P7ad 0.64 • IQ-'^ rad 0.78 ■ IQ-'* rads"' 1.16- 10-'" rads"' 0.92- lO-'-'rads-' 

ERM exp 0.58 • 10"' rad 0.69 • 10"' rad 0.52 • lO"' rad 0.77 ■ lO^'" rads"' 1.15 • IQ-'" rads"' 0.81 ■ lO"'-* rads"' 

gsf2006c 0.55 ■ 10-" rad 0.57 ■ lO'" rad 0.42 ■ lO''' rad 1.89 ■ lO''" rad s'' 2.00 ■ lO'" rad s'' 1.52 ■ 10-'" rad s'' 
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Fig. 6. The time series of the estimates of the daily offsets of nu- 
tation in obliquity when the MHB2000 nutation expansion was 
used as the a priori. The wrms is 9.8 ■ lO"'" rad. 
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rad s 

Fig. 7. The power spectrum of the estimates of the quasi-diurnal 
variations of components 1 and 2 of the perturbational vector of 
the Earth's rotation from solution A in the vicinity of the fre- 
quency of the near-diurnal free wobble. The estimate for the fre- 
quency -7.312026 ■ 10"^ rad s"', which corresponds to the tidal 
frequency is not shown. 



tide, turned out rather broad. The excerpt of the power spectrum 
produced from estimates of sine and cosine amplitudes of the 
components 1 and 2 is shown in Fig. [8] This signal cannot be 
attributed to the spectral leakage, since no excessive power was 
found in the vicinity of even a stronger tide at the frequency. 
A relatively broad-band signal in the vicinity of the -3Q„ fre- 
quency, i.e. Kt„ was found at the 3rd component of the rotation 
vector. The excerpt of the power spectrum produced from esti- 
mates of sine and cosine amplitudes of the component 3 is shown 
in Fig.|9l A weaker signal in the estimates can also be revealed 
in the vicinity of the -4f2„ frequency. A similar signal can be 
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Fig. 8. The power spectrum of the estimates of the quasi- 
diurnal variations of components 1 and 2 of the perturba- 
tional vector of the Earth's rotation from solution C in the 
vicinity of the -K2 frequency. The estimate for the frequency 
-1.458423 ■ lO^* rad s"', which corresponds to the -Kt, is not 
shown. 

seen at prograde frequencies in the vicinity K2,K-i,,K\ at com- 
ponents 1 and 2Q. 

Another peculiarity of the spectrum are sharp peaks at fre- 
quencies +4/50,,, +6/5Q,, at a level of 2-7cr above the noise 
level. No convincing explanation was found, but it is suspected 
that this signal in the estimates may be an artifact caused by er- 
rors in modeling by analogy with a detection of a very strong 
signal in estimates of the harmonic constituents of the pertur- 
bational Earth's rotation from GPS time series at frequencies 
that are mul tiple to the diurnal fre quency: S \,S2,S^,S \, etc., 
reported by Rothacher et al.l (1200 1'). It should be noted that no 
non-tidal signal at ^3,54 frequencies is seen from analysis of 
VLBI group delays. 

Solution D did not reveal other missed harmonic signals in 
the diurnal band. 

5.3. Error analysis 

The formal uncertainties of the amplitudes on harmonics con- 
stituents can be evaluated on the basis of the signal-to-noise ra- 
tio of fringe phases used for computing group delay by invoking 
the law of error propagation. These uncertainties are in a range 
of 5-12 prad. Analysis of the estimates of the constituents at the 
frequency bands where no tidal or no-tidal signal was detected 
provides a more reliable measure of noise in adjustments. It is 
16 prad for components 1,2 and 13 prad for component 3 for the 
diurnal band; 13 prad and 10 prad for these components at other 
frequency bands. This corresponds to displacements of 0.06- 

" Since component 3 was considered as a real value process, its spec- 
tral power at negative frequencies is the same as at positive frequencies. 
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Fig. 9. The portion of the power spectrum of estimates of the ter- 
diurnal variations of components 1 and 2 of the perturbational 
vector of the Earth's rotation vector in the ter-diurnal band. The 
broad peak is seen near the -Kt, frequency. 

0.12 mm at the Earth's surface. Evaluation of the level of sys- 
tematic errors is more problematic. The major possible source 
of systematic errors is considered to be a residual motion of the 
individual stations. In fact, the rotation of the station polyhedron 
was evaluated, and it was assumed that the motion of this polyhe- 
dron is a representative measure of the Earth's rotation. This as- 
sumption is valid to the extent that residual horizontal motion of 
individual observing stations is neglig ible. iPetrov & Mai (l2003h 
estimated harmonic site position variations and found that the ac- 
curacy of modeling the horizontal motion of individual stations 
is at the level of 0.4 mm. In the case of the errors of modeling be- 
ing completely uncorrelated, this error will be diluted as ^Nefj, 

where ^jNeff is the effective number of observing stations, 10- 
44, depending on how to define the effective number of stations. 
Unfortunately, the distribution of residual motions of stations at 
tidal frequencies shows a pattern of a systematic behavior, which 
does not support the hypothesis of uncorrelated errors. A con- 
servative estimate of the possible contribution of the unmodeled 
residual motion of the network of stations to the estimates of har- 
monic constituents of the perturbational rotation suggests a di- 
lution factor of 2, i.e. the surface displacements ~ 0.2 mm. That 
means systematic errors may be two times greater than random 

errors. 

iDehant et al.l (^003) investigated the influence of systematic 
errors due to the neglect of the modeling source structure. It was 
suggested to split the observed radio sources into two classes, 
"stable" and "unstable", and either to remove unstable sources 
from analysis or to estimate the time series of their positions. 
In this paper a different approach was used: proper motion of 
those sources that had a long enough history of observations was 
estimated. This method is supposed to reduce systemic errors in 
the estimates of the harmonic constituents in the perturbational 
rotation vector. 

6. Conclusion 

It was demonstrated that the empirical Earth rotation model can 
be determined directly from observations over a period of 22 
years using the least square estimation technique. The advan- 
tage of the proposed approach is that a continuous function de- 
scribing the Earth's orientation is determined in one step with- 
out producing intermediate time series. The consistency between 
station positions, source coordinates, and the empirical Earth ro- 



tation model is automatically achieved. Another advantage of 
the proposed approach is that a simplified a priori model with 
only 3 1 numerical parameters is sufficient, while the traditional 
approach needs a complicated a priori model of precession, nuta- 
tion, high frequency harmonic variations of the Earth's rotation, 
and a filtered and smoothed time series of the Earth orientation 
parameters produce d in the previous analysis, in total 46 000 nu- 
merical parameters jMcCarthv & Petit. eds.ll2064 ). 

The traditional approach to describing the Earth's rotation 
follows the f orm alism of either Newcomb and Andoyer or 
Guinoi d 19791) andlCapitaine et alj (|T986), and involves such no- 
tions as the celestial intermediate pole, the point of the vernal 
equinox, the non-rotating origin, the ecliptic, and other axes, 
points, planes, and circles on the celestial sphere. The advantage 
of the empirical Earth rotation model is that it is conceptually 
simpler, since it is built entirely kinematically and does not re- 
quire introduction of intermediate points, axes, planes that are 
not observable. 

It was demonstrated that the empirical Earth rotation model 
derived from analysis of VLBI observations gives the differences 
with respect to the EOF derived from analysis of independent 
GPS observations at moments of observation at the same level, 
within 20%, as the differences of the VLBI EOP series produced 
with the traditional approach. The advantage of the proposed 
approach is that the estimates of the EOP rates are a factor of 
1.5-2.0 closer to the GPS time series than the VLBI EOP rates 
estimated following the traditional approach. 

When results of analysis of observations are compared with 
theoretical predictions, two approaches can be taken: a) the pa- 
rameters that describe empirical data are formulated through pa- 
rameters of the theoretical models; b) theoretical predictions are 
transformed to a form that can be unambiguously determined 
from the observations. Representation of the Earth's rotation in 
the form of the expansion into basis functions establishes a foun- 
dation for the second approach. 

Scientific interpretation of the results of estimation of the 
empirical Earth rotation model will be given in the next paper. 

Acknowledgements. This work was done while the author worked at the 
National Observatory of Japan, Mizusawa, as a visiting scientist. The author 
would like to thank S. Manabe, T. Sato, Y. Tamura, D. Rowlands, Ch. Bizouard 
and an anonymous referee for useful discussion that helped to improve the 
manuscript. 

References 

Agnew, D. C, 1997, J. Geophys. Res., 102, 5109 

Aoki, S., Kinoshita, H.; Guinot, B., 1982, a, 105, 359 

Bizouard C, Brzezinski, A., Petrov S., 1998, J. of Geodesy, 72, 561 

Capitaine, N., Guinot, B., & Souchay, J., 1986, Celest. Mech., 39, 283 

Dehant, V., 1993, Adv Space Res., 13, 235 

Dehant, V., Feissel- Vernier M., de Viron O. et al., 2003, J. Geophys. Res., 108, 
2275 

Dickman, S. R, 1993, Geophys. J. Int., 112, 448 

Dick W. R. and B. Richter ed., 2005, lERS Annual Report 2004 of the 

International Earth Rotation and Reference Systems Service, Frankfurt am 

Main: Verlag des BKG, 152 p. 
Eubanks, T. M., 1993, Contributions of Space Geodesy to Geodynamics: Earth 

Dynamics: Geodynamics Series, vol. 24, ed. by D. E. Smith & D. L. Turcotte, 

AGU 

Getino J. & Femndiz J. M., 2001, MNRAS, 322, 785 

Guinot, B. 1979, In Time and the Earth's Rotation, ed. McCarthy D. D., & 

Pilkington, J. D. (D. Reidel Publishing Company), 7 
Herring, T. A., Gwinn, C. R. & Shapiro, 1. 1., 1986, J. Geophys. Res., 91, 4745 
Kalnay E., M. Kanamitsu, R. Kistler, et al., 1996, Bullet. Amer Meteorol. Soc, 

77, 437 

Kopeikin, S. M. and Schafer, G., 1999, Phys. Rev D, 124002 

Kotelnikov, V. A., 1933, In Proceedings Materialy k I vsesouznomu s'ezdu po 

voprosam technicheskoj rekonstruktcii dela svyazi i razvitija slabotochnoj 

promyshlennosti 



12 



L. Petrov: The empirical Earth rotation model from VLBI observations 



Mathews, P. M., Buffett, B. A., Herring, T. A. & Shapiro, I. I., 1991, 

J. Geophys. Res., 96, 8219 
Mathews, P. M., 2001, Journal of the Geodetic Society of Japan, 47, 231 
Mathews P. M., Herring, T. A., & BuflFett B. A., 2002, J. Geophys. Res., 107, 

10.1203/2001JB1000390 
Matsumoto, K., T. Takanezawa, & M. Ooe, 2000, Journal of Oceanography, 56, 

567 

McCarthy D. & G. Petit, eds., lERS Conventions, 2004, Verlag des Bundesamtes 

fuer Katrografie und Geodaesie, 127 
Moritz H., Earth Rotation: Theory and Observation, 1987, Frederick Ungar, 617 
NieU, A & Petrov L., 2004, In Proceedings of the workshop 'The State of 

GPS Vertical Positioning Precision: Separation of Earth Processes by Space 

Geodesy', ed. by T. van Dam and O. Francis, Luxembourg, 55 
Niimberger, G., 1989, Approximation by spline functions. Springer- Verlag 
Petrov L., & C. Ma, 2003, J. Geophys. Res., 108, doi: 10.1029/2002JB001801, 

2190 

Petrov L., & Boy, J.-R 2004, J. Geophys. Res., 109, 10.1029/2003JB002500, 
B03405 

Petrov L., 2005, In Proceedings of the 17th working meeting on European VLBI 

for geodesy and astrometry, ed. by M. Vennenbusch & A. Nothnagel, 113 
Ray, R. D., 1999, NASA/TM- 1999-209478, 58 

Rothacher, M., G. Beutler, R. Weber, and J. Hefty, 2001, J. Geophys. Res., 106, 
13,711 

Saastamoinen, J., 1972, In The Use of Artificial SatelUtes for Geodesy, Geophys. 

Monogr., AGU, 15, 247 
Schonberg, 1. J., Contributions to the problem of approximation of equidistant 

data by analytic functions, 1946, Quart. Appl. Math., 4 
Seideknann, P. K., 1992, Explanatory Supplement to the Astronomical Almanac 
Simon, J. L. & Bretagnon, P., Chapront, J. & Chapront-Touze, M., 1994, A&A, 

282, 663 

Souchay J. & Kmoshita H., 1996, A&A, 312, 1017 
Souchay J. & Kmoshita H., 1997, A&A, 318, 639 

Souchay J., Loysel B., Kinoshita H., and Folgueira M., 1999, A&AS, 135, 111 
Sovers, O. J, J. L. Fanselow, and C. S. Jacobs, 1998, Review of Modem 
Physics,70, 1393 

International VLBI Service for Geodesy and Astrometry, 1999 annual report, ed. 

N. Vandenberg, NASA/TP-1999-209243, Greenbelt, 308 
Wahr, J. M., 1981, Geophys. J. R. Ast. Soc, 64, 705 

Yseboodt, M., de Viron O., Chin T. M., & Dehant V, 2002, J. Geophys. Res., 
107, 2036, doi:10.1029/2000JB000042 



L. Petrov: The empirical Earth rotation model from VLBI observations, Online Material p 1 



Online Material 



L. Petrov: The empirical Earth rotation model from VLBI observations. Online Material p 2 
Table 3. Numerical values of the a priori Earth rotation model parameters used in data reduction. 



Var 
foi 

C02 

000 

^01 

002 

ZO 

Zl 

Z2 

eoo 

eoi 
eo2 
So 

£ln 

Pi 

P2 

ei 

ttl 

Q'2 
Pi 
P2 

Eo 
El 

E2 



7i 

72 



Value 
1.140216587056520 
3.542805701761733 
1.471291601425477 
9.909515599113584 
3.079019263961936 

-2.076601527511399 
1.140216587060519 
3.542805701761733 
5.331975251279779 
0.409092629687089 

-7.191223191481661 

-7.399638794037328 
1.753368559233960 
7.292115146706979 

-8.377867467753367 

-6.193374542381407 
4.473817016047498 
2.682642812740089 
2.182438855728973 
3.506953516079786 

-1.069696206302000 
3.982127698995000 
2.260937669429621 
1.029854567486117 

-7.875297448491237 
9.776692309499138 

-6.857935725000193 
3.783804480256964 
2.878954568890594 

-1.069696206302000 

-1.183000000000000 
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"Tir 



10-12 
10-" 

io-'2 

10-2^ 
lO-'o 

io-'2 

10-2^ 

10-'^ 
10-2' 

10-5 
10-5 
10-^ 
10-5 
10-^ 



io-« 

10-^ 

10-3 

io-'2 

10-22 
10-5 

10-** 
10-'^ 
10-^ 
10-'* 
10-** 



Units 
rad 

rad s- 1 
rad s-2 
rad 

rad S-' 
rad s-2 
rad 
rad S-' 
rad s-2 
rad 
rad S-' 
rad s-2 
rad 

rad S-' 

rad 

rad 

rad 

rad 

rad 

rad 

rad S-' 
rad S-' 
rad 
rad S-' 
rad s-2 
rad 
rad 
rad 
rad 
rad S-' 
rad S-' 



(1994) 
(1994) 
(19941 
. (1994) 
(1994) 
. (1994) 
, (1994) 



Source 

Simon etal. (1994) 
_Simonetal. (1994) 
Simon etal. (1994) 
Simon etal. (1994) 
Simon et al. (1994) 
^imon et al^ 
Simon et al.^ 
Simon et al.^ 
Simon et al. 
Simon et al. 
Simon et al. 
Simon et al. 
Aoki (1982) 
Aoki ( 1982) 

Souc hav & Kinoshita (1996) 
Souchav & Kinoshita (199&) 
Souchav & Kinoshita (1996^ 
Souchav & Kinoshita (1996) 
Simon et al. (1994) 
Simon etal. (1994) 
Simon et al. (1994) 
Simon et al. (1994) 
LSQ fit 
LSQfit 
LSQ fit 
Dickman 
Dickman 
LSQfit 
LSQfit 
Dickman 
LSQ fit 



